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Abstract 

Gaussian fields (GFs) are frequently used in spatial statistics for their versatility. 

The associated computational cost can be a bottleneck, especially in realistic appli¬ 
cations. It has been shown that computational efficiency can be gained by doing the 
computations using Gaussian Markov random fields (GMRFs) as the GFs can be seen 
as weak solutions to corresponding stochastic partial differential equations (SPDEs) 
using piecewise linear finite elements. We introduce a new class of representations 
of GFs with bivariate splines instead of finite elements. This allows an easier im¬ 
plementation of piecewise polynomial representations of various degrees. It leads to 
GMRFs that can be inferred efficiently and can be easily extended to non-stationary 
fields. The solutions approximated with higher order bivariate splines converge faster, 
hence the computational cost can be alleviated. Numerical simulations using both 
real and simulated data also demonstrate that our framework increases the flexibility 
and efficiency. 

Keywords: Gaussian Markov random field; Spatial approximation; Multivariate splines; 

Non-stationary spatial process; Mapping. 
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1 Introduction 


Gaussian fields (GFs) are at the core of spatial statistics, especially in the class of struc¬ 
tured additive regression models, named latent Gaussian models, which are flexible and 


extensively used (Gressie, 1993 Banerjee et ah, 2004; Diggle and Ribeiro, 2007). However, 


when making statistical inference, it is usually needed to evaluate the likelihood function 
or the latent Gaussian held distribution, for which we need to make computations on dense 
matrices, e.g. the covariance matrix S(0), typically of order 0{'n?) where n is the dimen¬ 
sion of S(0). Rue et ah (2009) overcome this computational hurdle. They approximate 


Bayesian inference in latent Gaussian models by assuming that the latent held is Gaus¬ 
sian Markov random held (GMRF). With only a few hyperparameters, integrated nested 
Laplace approximations (INLA) produce faster inference than simulation based approaches 


such as MGMG. To take advantage of the computational efficiency of GMRF, Lindgren et ah 


(2011) constructed an explicit link between GFs and GMRFs. They considered the GFs 


with Matern covariance function. 
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r u,v = 


2^-iF(i/) 


K V 
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K V 
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( 1 ) 


where ||v — u|| is the Euclidean distance between two locations u and v G is 

the modihed Bessel function of the second kind and order u > 0, k > 0 controls the 
nominal correlation range through p = corresponding to correlations near 0.1 at 

the Euclidean distance p, and is the marginal variance. The integer value of u determines 
the mean-square diherentiability of the underlying process. They noticed that a Gaussian 
held a:(u) with Matern covariance ([^ is a solution to the linear fractional SPDE 


— A)"/^(ra;(u)) = hF(u), 


u G 




a = u + d/2, K > 0, z/ > 0, 


( 2 ) 


where the i nnovation process W is spatial Gaussian white noise with unit variance (Whittle] 


1954 


1963|), A = ^ is the Laplacian operator, and r controls the marginal variance 


through the relationship 


r2 = 


T{u) 


F(z/ -|- d/2){ATTY/'^K,'^^ 




Therefore to hnd the GF a;(u) with covariance function ([^ is to hnd the solution to 
1^. Denote the inner product of two functions / and g on as (/, g) = f{u)g{u) d u, 
then the weak stochastic solution to SPDE ([^ can be found by requiring that 


(k" - Ar'^Tx) i {,p, W), 


(3) 


for suitable functions 0(u), where ‘ = ’ denotes equality in distribution (Walsh, 1986). Then 


Lindgren et ah 

|2011) constructed a hnite element representation ( 

Brenner and Scott 

2008 


of the Gaussian random held over an unstructured triangulation of the form 


k=l 


( 4 ) 
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where {V’A:}fc=i are piecewise linear basis functions. They showed that the Gaussian weights 
{wfc}fc=i are GMRFs when a = 1, and can be approximated with GMRFs when a > 2. 
Therefore the computations for GFs can be carried out using GMRFs and the computa¬ 
tional efficiency can be improved dramatically. This work is closely related to the spatial 
spline regression models by Sangalli et al. (2013) where a spatial surface is approximated 
with hnite elements. Another recent related work is Nychka et al. (2014), where the authors 
proposed a representation of a random held using multi-resolution radial basis functions on 
a regular grid. They also assumed that the coefficients associated with the basis functions 
to be distributed according to a GMRF to speed up the computation. 

It is stated in Lindgren et al. (2011) and Simpson et al. (2012) that the convergence 
rate of a hnite element approximation to the full solution to the SPDE (|^ is of order 
0{h’^) where h is the length of longest edge in the triangulation. Hence the convergence 
can be achieved by rehnement of the underlying triangulation which is usually called the 
h-version hnite elements. An alternative is to increase the approximation order over any 
hxed triangulation with higher degree polynomials over each triangle, which is called the 
p-version hnite elements (the degree of polynomials is usually denoted by p). It has been 
illustrated that the convergence rate of the p-version cannot be worse than the h-version in 
most cases (Babuska et ah, 1981). To increase the approximation order over each triangle, 
multivariate splines over triangulations can be employed instead of conventional hnite ele¬ 
ments. This provides a hexible and easy construction of splines with piecewise polynomials 
of various degrees and smoothness. Basic concepts and theories of multivariate splines can 
be found in the monograph by Lai and Schumaker (2007). Multivariate splines have been 
shown to be more efficient and hexible than conventional hnite element method in data 


htting problems and solving PDEs, see 

Awanou et al. 

(2005 

). It has been applied in spatial 

statistics. For example. 

Guillas and Lai 

(2010 

) introduced a spatial data analysis model 


with bivariate splines by penalizing the roughness with a partial diherential operator; this 
has been demonstrated to be more efficient and accurate than thin-place splines in the 
application of ozone concentration forecasting (Ettinger et ah, 2012). In this paper, we 


introduce bivariate splines to represent the GFs on and show its advantages over the 
piecewise linear hnite elements used in Lindgren et al. (2011). Within our framework of 
the SPDE approach using bivariate splines, it is allowed to choose piecewise polynomial 
representations of arbitrary degrees to adapt to the various data structures and features, 
and make the inference more computationally efficient. 

The paper is structured as follows. In Section some basics of bivariate splines in 
the Bernstein form (B-form) are reviewed hrst. Then we show how to link the GFs with 
GMRFs within the framework of bivariate splines, establish the theoretical properties of the 
bivariate spline approximations and discuss extensions to non-st at ionary helds. In Section 
1^ we conduct several numerical simulations to illustrate our method and compare with the 
approach of Lindgren et al. (2011) on both real and simulated data sets. Sectionconsists 
of conclusion and discussion. Proofs are in the Appendix. 
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2 SPDE approach using bivariate splines 


2.1 B-form bivariate splines 

Let A be a triangulation of a bounded domain C We consider the continuous spline 
spaces 

5°(A) = {s e e iP,,VT e A}, 

where Vd is the space of bivariate polynomials of degree d > 1 , (7°(hi) is the space of all 
continuous functions on hi. For any d > 1, the spline space ^^(A) contains all possible 
continuous spline functions which are bivariate polynomials of degree d over each triangle 
T G A. We apply the B-form representation of splines in 5° (A) proposed by Awanou 


et ah (2005) in this paper, which allows an easy construction of locally supported basis 


functions and associated calculations. In this section we only give a brief introduction to 


the bivariate splines and some details are relegated to Appendix A.l For more complete 


and in-depth explanations, see Lai and Schumaker (2007). 


Let T = (vi, V 2 , V 3 ) be a non-degenerate (i.e. with non-zero area) triangle with vertices 
vi = (xi,?/i), V 2 = {x 2 ,y 2 ) and V 3 = ( 0 : 3 , 2 / 3 ). Then every point v = {x,y) G has a 
unique representation in the form 


V = 61V1 -h b2'V2 hvs, 


(5) 


with bi + b 2 + bs = 1 , where 61 , &2 and 63 are named the barycentric coordinates of the 
point V = (x, y) relative to the triangle T. The polynomials 




U bP b^ 


( 6 ) 


are called the Bernstein polynomials of degree d relative to triangle T. Then for each spline 
function s G S'°(A), we can write 


i>iT= 5^ cS.B.yt Tea, 

2+J + fc = cZ 


where the coefficients c = {c7^, i + j + k = d,T G A} are called B-coefficients of s. 


2.2 SPDE modelling with B-form bivariate splines 

Let {' 01 , V’ 2 , be a set of locally supported basis functions of Sj^(A) for any d > 1, 

where m = dimS'°(A) (see Appendix A.l and Lai and Schumaker (2007) for more details 


on these basis functions). For any h = 1, the corresponding B-coefficients of ifjh are 

denoted by c^. Then the weak solution to the SPDE (|^ in the spline space S^(A) can be 
represented as 


XaIu)= 




(7) 


h=l 


Remember that the weak solution to SPDE ([^ can be found by requiring ( [^ for any 
sensible test functions. However it is not possible to test all functions. Following Lindgren 




















et al. (2011) we can choose a finite set of test functions. Specifically, we choose (t)h = 
— A)h^/i for a = 1 leading to the least squares solution. For a = 2, we can choose 
either (ph = i^h for any d > 1 or (ph = — A)'0/i for d > 2, leading to the Galerkin or 

least squares solution respectively. For a > 3, if we let a = 2 on the left-hand side of 
the SPDE (I 2 I) then the right-hand side is a Gaussian process generated by the operator 


Then we can choose (ph = i^h for this innovative SPDE. Hence we get a 


recursive Galerkin solutions ending with either a = 
results as below. 


1 or 2. Therefore we have the main 


Theorem 1 The vector of weights w = {wh, h = 1, ...,my of bivariate spline representation 
for the solution to SPDE ([^ defined in ([^ is Gaussian with mean zero and the precision 
matrix Qq are given as follows: 

(1) for a = 1, 

Qi = r2(/s:2M + K), 

(2) for a = 2, 

qG = + 2k2k + KM-^K), 

+ 2k^¥. + R), 

where and are the Galerkin and least sguares solutions respectively, 

(3) for a > 3, 


Qa = «"Qa-2 + k 2(Q«-2M-1K + KM-^Q„_2) + KM-^Q^.sM-^K, 


where 

M = C'MoC, K = C'KoC, R = C'RoC, 

and Mo = diag(Mr,T G A), Kq = diag(KT,T G A) and Ro = diag(RT,T G A) are 
block diagonal sguare matrix with sguare blocks 


IVIt’ — 



B'[jkix,y)B^^^{x,y)dxdy 


U^fL-\-K=d 


i-\-j+k=d 


and 


Kt = 


R 7 


2 /) y)dxdy 


^B'kki^, y)AB^^^{x, y)dxdy 


IT 


h'+/i+R=d 

i-\-j+k=d 

l'-\-fL-\-K=d 

i-\-j+k=d 


respectively and C = (ci,...,Cm) whose h-th column is the B-coefficient vector of basis 
function fjh- 
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Since the basis functions {-^i, ^’ 2 , •••, V’m} are locally supported in S'^(A), the matrices 
M, K and R are guaranteed to be sparse. However in the Galerkin solution for a = 2 
and recursive solutions for a > 3, the inverse matrix can be not sparse, making the 
precision matrix dense. The mass lumping technique (Chen and Thomee, 1985) can be 


applied by replacing M with a diagonal matrix M whose elements are the sum of each 
row of M, i.e. Mjj = 'Yhj Mjj. Therefore the precision matrix is sparse and the underlying 
coefficients w are approximated with a GMRF. In the next sections we show the effect 
of such Markov approximation in both theory and application. Another point we want 
to mention is that piecewise linear hnite elements are typical spline functions. The hnite 


element representation in Lindgren et ah (2011) is just a bivariate spline in *S'°(A 


2.3 Approximation properties 

Previously, we have constructed the bivariate spline approximation a;A(u) to the true Gaus¬ 
sian random held or the solution to the SPDE ([^. We explore some theoretical properties 
of our method in this section. 

Dehne the Hilbert space associated with the differential operator {k^ — A) to be 
the space of square integrable functions f{x,y) for which 


V f{x, y) ■ V f{x, y)dxdy is hnite following 

Lindgren et al. 

(2011 

). Approximation results 

for bivariate splines, e.g. Th. 5.19 in 

Lai and Schumaker I 

2007), show that the bivariate 


spline space S'^(A) for any d > 1 spanned by a hnite set of basis functions {'ipi, ...,■0^} 
is dense in H^: for every / G H^, there is a sequence {fm}, fm £ >S'°(A) such that 
limm->.oo \\f — fm\\m =0 where the limit scenario m —)■ cx) corresponds to |A| —)■ 0 where 
|A| is the length of the longest edge in the triangulation A. Using this fact, it follows 


directly from the Th. 3-4 in Appendix G.2 of Lindgren et ah (2011) that, the bivariate 


spline approximation xa converges weakly to the weak solution to the SPDE. Note that 
the weak convergence of a:A obtained for cannot be derived directly but can be easily 
proved in the same fashion with just a few modihcations. In addition, we can derive rates 
of convergence results. For example, when a = 2 we have the proposition below regarding 
to the Galerkin solutions. 


Proposition 1 Let L = — A), x/s.{s) is the bivariate spline approximation of the 

random Gaussian field x{s) in the spline space S'°(A), d > 1. Then for any f G f l 
with 1 < m < d, where is a Sobolev space that is defined in Appendix\A.t^ 


we have 


E 


f{s)L{x{s) - x^{s))ds ] < iF|Ar+i|/U+i, 2 ,o 


where K is a constant, |A| is the length of the longest triangle edge in the triangulation A 


and |/|m+i, 2 ,o is defined in Appendix A.3 


It is clear that we are able to achieve a faster convergence rate by using bivariate splines 
with higher degree d. For example, when d = 3 the convergence rate can be as high as 
(P(|A|^), which is two magnitude higher than (P(|Ap) 


m 


Lindgren et ah (2011). 


As we have mentioned, the matrix M in the Galerkin solutions is lumped by replacing 
M with a diagonal matrix M which yields a Markov approximation £ a to the bivariate 
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spline solution xa- Let / and g be test functions in and let /a and ^^a be their 
projections onto the bivariate spline space S^{A) for any d > 1, with basis weights w/ 
and Wg. Since the recursive algorithm for a > 3 is based on a = 2 at each iteration, here 
we only investigate the effect of the Markov approximation on the Galerkin solutions for 
a = 2. When a = 2, the difference between the covariances for the Markov approximation 
xa and the bivariate spline solution xa is 

eA(/A, 5'a) = cov((/, Lxa)o, {g, Lxa)o)-cov((/, Lx^)n, {g, Lx^)n) = w'j-Mwg-w'j-Mwg. 
We have the following result showing that such a difference can be bounded. 

Proposition 2 For /a, ^A e Sj^(A), we have 

|eA(/A,5'A)| < iL|Ap 


where K is a positive con stant dependent on ||/A|| 2 , 2 ,n, Us'Alhg.n, il/A||oo,o, ||fi'A||oo,n (norms 
are defined in Appendix A.3} and |A| is the length of the longest triangle edge in the 
triangulation A. 


We can see that the convergence rate of the bivariate spline representation in *S'°(A) in 
Proposition!^ may be decreased. But it is still at least (P(|Ap) in theory, which is as good 
as the linear hnite elements representation. In fact, numerical simulations later illustrate 
that the approximation using bivariate splines with higher degree d can be more efficient. 


2.4 Non-stationary fields 

showed that the SPDE (|^ can be extended to a non-stationary 

(k=(u) - A)“/^(r(u)x(u)) = W{u), (8) 

where the parameters and r are not constants but depend on the location u. The two 
parameters are assumed to vary slowly over the domain of u and have the general form of 
low dimensional representation 


Lindgren et al. (2011 


version 


\og{K\u)) = 

1=1 




U 


log(r(u)) = 

1=1 


where the number of smooth basis functions { B ^ ^(u)} for k^(u) and t(u) respectively, 
and Ur, should not be large to guarantee computational efficiency. The inner product can 
be approximated with 

{'fit, ~ l'i‘^{u*){'fit,'fis), 

where u* is some point in the support of fig which can be chosen to be the domain point 
associated with the non-zero B-coefficients of fig (see Appendix |A.1| for more details about 
domain points and their relationship to fig). Defining the diagonal matrices 

Kf = diag(K^(^fe), h = 1,..., m), r = diag(r(^fe), h = 1,..., m). 
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where is the domain point associated with the non-zero B-coefficients of basis function 
'iph for h = 1, It can be easily shown with minor modihcation of the proof of Theorem 

[^that the weights w in the bivariate spline representation Q can be approximated with 
GMRF as well. For example when a = 2, the precision matrix of w is 


t) = + 2k,^K + KM^^K)t, 

0,2^t) = -I- 2k^K + R)t, 


for Galerkin and least squares solutions respectively. As stated in Lindgren et ah (2011), 
by assuming the parameters and r to be constant locally, the solution to the SPDE 


can still be interpreted as a Matern held over a local area and the associated global 
non-stationary held can be achieved by combining all the local Matern helds automatically 
via the SPDE. 


3 Numerical simulations 


We conduct several numerical simulations to evaluate the performance of the SPDE ap¬ 


proach with bivariate splines and compare with the linear hnite element approach in Lind- 
gren et ah (2011) in terms of spatial prediction. In all simulations over we hx a = 2 


which corresponds to the smoothness parameter z/ = 1 in the Matern covariance function. 
The full Bayesian inference for the model is run in R-inla (www.r-inla.org) using the 


integrated nested Laplace approximation (Rue et ah, 2009). For brevity, our proposed bi 


variate spline approximation in *S'°(A) is denoted BS-SPDE with d = 1, 2,... (BS-SPDE-G 
or BS-SPDE-LS for Galerkin or least squares solution respectively) while the linear hnite 
element approximation is denoted LFE-SPDE. 


3.1 Comparison of LFE-SPDE and BS-SPDE 

Study 1 

In this simulation, we compare the LFE-SPDE method and BS-SPDE-G or BS-SPDE- 
LS of degree d > 2 in data htting and prediction for some common surfaces. Elevations 
of diherent surfaces are collected on a grid over square [—2,2] x [—2,2] that is equally 
spaced every 0.2. We consider the specihc surface to be an unknown random Gaussian 
held with Matern covariance which is also the solution to the stationary SPDE (|^. The 
whole random held can be approximated using LFE-SPDE and Galerkin or least squares 
BS-SPDE respectively. Then elevations on another hner grid that is equally spaced every 
0.01 over square [—2, 2] x [—2, 2] can be predicted. The prediction accuracy for the whole 
surface can be measured with mean squared error MSE = X]r=i(/(w) — /(w))^/''^, where 
/(u*) is the true elevation on location u* and /(uj) is the prediction using corresponding 
posterior means. 

We consider four diherent surfaces here including 2sin(x) cos{y) and 2exp(—) with 
three diherent shape parameters s = 2, 1, 0.5. As we have shown that both the hnite 
element and bivariate spline approximations converge weakly to the full SPDE solution. 
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we construct 35 different meshes that have 2, 3, 6, 13, 19, 28, 53, 96, 112, 148, 212, 279, 
342, 390, 444, 520, 705, 874, 1065, 1368, 1802, 2416, 2798, 3176, 3708, 4428, 5514, 6696, 
8460, 10958, 15009, 21832, 26718, 33776, 43875 triangles respectively to demonstrate the 
convergence (mesh size |A| monotonically decreases roughly from 6.6 to 0.026). For each 
approach, the associated number of basis functions (denoted by Nb) and CPU time for 
inla (denoted by T^pu in seconds) are recorded when the corresponding MSEs reach levels 
of 10“^, 10“^, 10“^, 10“^, 10“^, 10“®, 10“^, 10“®. Nb is also the dimension of corresponding 
precision matrix of the weights w and directly relates to the computational complexity in 
the associated calculations. For example the samples and likelihoods can be computed in 
operations for two dimensional GMRFs. For comparison, the simulation stops 
when the number of basis functions of BS-SPDE with d > 2 exceeds the number of basis 
functions of LFE-SPDE using the densest mesh. The results are presented in Figure 
where the y-axes for Nb and Tcpu are taken on a logarithmic scale. 

From Figure [T] we can see that in general BS-SPDE with d > 2 can be more efficient 
than LFE-SPDE both in terms of number of basis functions and computing time needed 
to reach specihc levels of MSE, especially those lower than 10“^. Note that in the left 
side of Figure [T| the dash lines for BS-SPDE-LS are invisible as they coincide with the 
solid lines for BS-SPDE-G while BS-SPDE-LS is more computationally efficient in general 
than BS-SPDE-G as shown in the right side when the associated precision matrices are of 
the same dimension. Specihcally, for the surface 2 sin(a;) cos(|/), to reach a specihc level 
of MSE, especially high precision levels such as 10“® or 10“^, BS-SPDE-G and BS-SPDE- 
LS with high degree d > 3 are more efficient since they require only less than 10% of 
the basis functions and computing time required by LFE-SPDE. For the Gaussian surface 
2exp(—), BS-SPDE with d>2 are generally much more efficient than LFE-SPDE for 
the MSE levels up to 10“^ with about 50% gains in the computing time. But BS-SPDE-LS 
does not reach the MSE level 10“® and BS-SPDE-G with d = 3 takes more computing 
time than the others to reach the MSE level 10“®. For the next Gaussian shape surface 
2exp(— ^ ) which is steeper than the previous one, BS-SPDE with high degrees can 

be better than LFE-SPDE for the MSE levels around 10“^ to 10“® but their efficiency 
is decreased to reach the higher precision levels 10“^ and 10“®. Only BS-SPDE-G with 
d = 2 reaches the lowest MSE level 10“®. However, BS-SPDE-LS with <7 = 4 reaches the 
low MSE level 10“^ within only 20% of the computing time required by LEF-SPDE. For 
the last surface which is quite steep, even though neither BS-SPDE-G nor BS-SPDE-LS 
with d > 2 is more efficient than LFE-SPDE in most cases, BS-SPDE-G with d = 2 is 
comparable with LFE-SPDE and reaches the high precision levels 10“®, 10“^ by requiring 
slightly less number of basis functions and similar time, and BS-SPDE-LS with d = 4 is 
more efficient than the others to reach the MSE level 10“®. 

From these results, we can conclude that BS-SPDE can be much more efficient in many 
cases especially when the high precision levels are desired and the target functions are 
smooth. For functions that are not that smooth, lower degree representations such as 
LFE-SPDE or BS-SPDE with d = 2 might be more appropriate, which is also consistent 
with the general comments by Babuska et ah (1981). Note that even for the last Gaussian 
shape surface which is much less smooth than the others, BS-SPDE-G with d = 2 still can 
be comparable with LFE-SPDE; and we still obtain 50% gains in the computing time using 
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2sin(x)cos(y) 



LFE-SPDE 
BS-SPDE-Gd=2-B ■ BS-SPDE-LS d=2 
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BS-SPDE-Gd=4-A- BS-SPDE-LS d=4 
BS-SPDE-G d=5 BS-SPDE-LS d=5 a 
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2exp((x^ + y^)/l) 
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2exp((x^ + y^)/0.5) 



MSE MSE 


Figure 1: Number of basis functions {Nb) and CPU time for inla {Tcpu in seconds) required 
by LFE-SPDE, BS-SPDE-G and BS-SPDE-LS with d = 2,3,4, 5 respectively to reach 
specihc MSE levels for different surfaces 
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BS-SPDE-LS with d = 4 if the MSE level 10 ® is desired. 


Study 2 


In this study, we compare LFE-SPDE and BS-SPDE in spatial estimation and prediction 


with real data sets that are extracted from the ETOPOl Global Relief Model (Amante and 


Eakins, 2009), which is a 1 arc-minute global relief model of Earth’s surface that integrates 


land topography and ocean bathymetry. The data is available from National Geophysical 
Data Genter (NGDG), USA. Four different regions around the the Strait of Juan de Fuca 
area are chosen for this study as shown in Figure In general, region 1 covers relatively 
simple and gradual variations in the near shore seabed while the seabed in the other three 
regions is quite complicated. 



Figure 2: Four regions extracted from ETOPO 1 around the Strait of Juan de Fuca from 
NGDG 


For comparison, both in-sample and out-of-sample predictive £t performance are ex¬ 
plored using LFE-SPDE, BS-SPDE-G and BS-SPDE-LS with d = 2,3,4 based on various 
meshes. We denote the observations by |/i, 1 / 2 , ..., |/n- As for the in-sample £t measure¬ 
ment, root mean square error (RMSE) between the observations and the predictions at the 
observed locations 

n 


RMSE 




i=l 
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are calculated where the predictions jji are taken to be the associated posterior mean. 
Since the SPDE approach aims to estimate the whole surface, smaller RMSE suggests the 
estimated surface is closer to the measurements at the observed locations. To measure 
the predictive performance, leave-one-out cross validation is employed using the embedded 
function within R-inla. The logarithmic score (Log Score) of prediction is dehned as 


N 


Log Score = - — ^ log[7r(?/i|?/_*)], 


i=l 


where 7r(?/j|?/_j) is the posterior predictive density of Ui given all the other observations 
U-i. Therefore the smaller Log Score is, the more certain we are with the predictions. 
Furthermore six meshes are built which are denoted by mesh 1-6 respectively as shown 
in Figure The meshes are extended with coarse triangles to avoid boundary effect 
(Lindgren et al., 2011). The number of basis functions for each combination of mesh and 


SPDE method is shown in Table[^ Table [^presents the RMSEs and Log Scores using LFE- 
SPDE, BS-SPDE-G and BS-SPDE-LS with d = 2,3,4 based on the six meshes respectively. 



mesh 4 



mesh 2 



mesh 5 



mesh 3 



mesh 6 



Figure 3: Six meshes for Study 2 

In general, it can be shown from Table that as the triangulation becomes denser, 
the estimations and predictions are more accurate using both LFE-SPDE and BS-SPDE- 
G/LS in most cases. For a particular mesh, the RMSEs and Log Scores of BS-SPDE-G/LS 
with d > 2 are generally smaller than those of LFE-SPDE. In terms of the number of 
basis functions, BS-SPDE-G/LS with d > 2 also demonstrate better performance than 
LFE-SPDE in most cases. For example, for region 4, the RMSEs and Log Scores of BS- 
SPDE-G with d = 2 based on mesh 4 or BS-SPDE-G with d = 4 based on mesh 2 are 
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Table 1: Number of basis functions using LFE-SPDE and BS-SPDE-G/LS with d = 2,3,4 
based on six meshes 


mesh 

LFE-SPDE 

BS-SPDE d = 2 

BS-SPDE d = 3 

BS-SPDE d = 4 

1 

28 

89 

184 

313 

2 

70 

252 

547 

955 

3 

195 

749 

1663 

2937 

4 

244 

945 

2104 

3721 

5 

536 

2111 

4726 

8381 

6 

978 

3881 

8710 

15465 


much smaller than corresponding components in LFE-SPDE based on mesh 6 while they 
have similar numbers of basis functions. BS-SPDE-LS approach shows similar properties 
with BS-SPDE-G in general. We notice that BS-SPDE-LS yields smaller RMSEs than 
BS-SPDE-G in most cases. However, in terms of Log Score, BS-SPDE-G performs better 
than BS-SPDE-LS in most cases for the other three regions except region 1. Another point 
we want to mention is the sudden change in the model performance due to the rehnement 
of mesh or increase of polynomial degree d. For example, the RMSE obtained using LFE- 
SPDE for region 2 is about 15.87 based on mesh 4; it is decreased suddenly to only 0.037 
based on mesh 5 or 0.012 using BS-SPDE-G with d = 2. We notice that the associated 
number of basis functions increases. This sudden change may because the hnite elements 
or splines reach some level of degree of freedom that is enough to model the surface well. 

Based on Table we can also select the models with good performance for continuous 
map reconstruction among the different combinations of meshes and SPDE approaches for 
each data set. Note that in most cases it is difficult to have a model with smallest RMSE 
and Log Score at the same time, so we only choose the one with relatively small RMSE and 
Log Score. In this way, the reconstructed map can be close to the elevations at the observed 
locations; meanwhile we are more confident with the predictions at the other locations. As 
marked with asterisks in Table we choose BS-SPDE-LS with d = 4 based on mesh 3 for 
region 1, BS-SPDE-LS with d = 3 based on mesh 3 for region 2, BS-SPDE-G with d = 4 
based on mesh 3 for region 3, and BS-SPDE-G with <7 = 4 based on mesh 3 for region 
4. Note that for region 1, BS-SPDE-LS with d = 4 based on mesh 6 yields both smallest 
RMSE and Log Score among all the combinations. However the associated computational 
cost is much heavier than the others. There is some trade off between model performance 
and computational cost. Hence we select the one with relatively good performance and is 
also computationally efficient. Then the posterior means and standard deviations of the 
four regions predicted using the respective selected models are displayed in Figure The 
posterior means in general capture the main features of the corresponding regions and the 
posterior standard deviations provide uncertainty estimates of the predictions. Note that 
the selection rule of predictive model here is quite simple and subjective. More appropriate 
model selection techniques can be employed in application. 
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Table 2: Study 2: RMSE and Log Score (RMSE|Log Score) using LFE-SPDE, BS-SPDE-G 
and BS-SPDE-LS with d = 2,3,4. The values in red colors are obtained based on similar 
number of basis functions between 945 to 978. The values marked with asterisks (*) are 
the selected model fit for map reconstruction 


mesh 

LFE-SPDE 

BS-SPDE-G 

BS-SPDE-LS 

d = 2 

d = 3 

d = 4 

d = 2 

d = 3 

d = A 

region 1 

1 

1.85|2.07 

1.41|1.84 

1.18|1.73 

1.11|1.73 

1.41|1.84 

1.20|1.72 

1.16|1.71 

2 

1.17|1.71 

0.76|1.53 

0.55|1.47 

0.29|1.21 

0.88|1.65 

0.81|1.70 

0.74|1.73 

3 

0.83|2.01 

0.46|1.48 

0.33|1.42 

0.16|1.09 

0.43|1.64 

0.048|1.18 

0.027|0.98* 

4 

0.79|1.61 

0.46|1.48 

0.34|1.43 

0.17|1.16 

0.23|1.43 

0.039|1.08 

0.028|0.99 

5 

0.62|1.51 

0.49|1.50 

0.42|1.45 

0.36|1.38 

0.024|1.51 

0.022|1.30 

0.021|1.08 

6 

0.63|1.52 

0.50|1.51 

0.49|1.48 

0.44|1.53 

0.021|1.74 

0.020|1.19 

0.019|0.92 

region 2 

1 

79.81|5.82 

66.01|5.66 

52.26|5.41 

47.01|5.30 

66.41|5.65 

52.21|5.42 

45.81|5.27 

2 

36.20|5.02 

21.27|4.53 

9.09|3.86 

0.02|0.50 

22.47|4.58 

14.52|4.23 

12.48|4.10 

3 

17.24|4.34 

0.013|0.51 

0.0091|0.43 

0.0088|0.45 

0.015|0.43 

0.011|0.41* 

0.0094|0.50 

4 

15.87|4.28 

0.012|0.52 

0.010|0.46 

0.010|0.50 

0.013|0.47 

0.0098|0.46 

0.0091|0.54 

5 

0.037|1.85 

0.013|0.85 

0.010|0.55 

0.010|0.48 

0.0086|0.61 

0.0069|0.55 

0.0065|0.57 

6 

0.032|1.98 

0.011|0.72 

0.010|0.55 

0.012|0.53 

0.0075|0.71 

0.0062|0.64 

0.0073|0.81 

region 3 

1 

41.88|5.18 

37.52|5.09 

30.83|4.89 

25.44|4.69 

37.52|5.10 

30.51|4.89 

25.21|4.68 

2 

26.64|4.72 

11.62|3.93 

6.06|3.43 

0.021|0.45 

11.70|3.94 

7.83|3.59 

0.029|0.52 

3 

11.45|3.95 

0.028|0.74 

0.016|0.50 

0.015|0.48* 

0.022|0.58 

0.021|0.56 

0.017|0.61 

4 

10.18|3.86 

0.025|0.72 

0.017|0.53 

0.015|0.48 

0.021|0.55 

0.019|0.60 

0.016|0.59 

5 

7.21|3.66 

0.021|1.01 

0.018|0.64 

0.016|0.50 

0.014|0.69 

0.012|0.68 

0.012|0.72 

6 

0.053|2.08 

0.019|0.99 

0.019|0.80 

0.019|0.64 

0.012|0.82 

0.011|0.79 

0.012|0.96 

region 4 

1 

47.12|5.30 

35.98|5.07 

28.46|4.83 

22.95|4.61 

36.02|5.07 

28.46|4.83 

22.73|4.59 

2 

26.36|4.71 

12.56|4.02 

8.08|3.71 

0.020|0.48 

12.47|4.01 

9.24|3.78 

0.033|0.55 

3 

13.85|4.13 

0.023|0.62 

0.016|0.52 

0.013|0.43* 

0.019|0.55 

0.016|0.54 

0.014|0.56 

4 

12.42|4.04 

0.021|0.76 

0.017|0.54 

0.014|0.46 

0.019|0.54 

0.015|0.56 

0.014|0.59 

5 

0.061|2.19 

0.019|0.93 

0.017|0.64 

0.015|0.49 

0.012|0.63 

0.011|0.63 

0.010|0.68 

6 

0.049|2.14 

0.017|0.93 

0.018|0.78 

0.017|0.58 

0.011|0.77 

0.010|0.77 

0.010|0.88 
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Figure 4: Posterior mean (top) and standard deviation (bottom) for the regions 1-4 (left 
to right), after selection of the appropriate approximation models (marked with asterisks 
in Table 

3.2 Spatial analysis of ozone levels data over Eastern USA 

In this section, we analyse a data set of ozone levels at a certain hour in one of days in 
September, 2005 around the Eastern United States, which is available from the Air Explorer 
Database of Environmental Protection Agency (EPA), using the non-stationary BS-SPDE- 
G method. The data set has 546 locations where ozone levels are recorded. As shown in 
Figure the observations of ozone concentration are distributed unevenly and the domain 
is irregular. Denote the ozone levels by Zi and the associated locations by Sj = 
i = I, ...,546. We consider a simple spatial model 

Zi^bo + f{si), i = l,...,546, 

where bo is the intercept and the spatial effect /(sj) is assumed to be a non-stationary GF 
generated by the non-stationary version SPDE ([^, approximated with bivariate splines in 
S'^(A) with d = 1, 2,3,4, 5. The triangulation A is shown in Figure 

The non-stationary parameters 'r(u) and k{u) are represented with two-dimensional B- 
splines that have Ux and Uy basis functions in the a:-direction and ^/-direction respectively. 
Therefore at any location s = {x,y), the basis functions of the associated B-spline can be 
calculated as Bik{s) = Bf{x)Bl{y), where i?f (•) and Bl{-) are the basis functions in x and 
y directions respectively, for I = l,...,nx and k = l,...,ny. Hence there are UxUy basis 
functions in total for each of the parameters k^{-) and r(-). We consider 12 models A-L 
with different combinations of the number of basis functions in Table [H Note that with one 
basis function, the B-spline is constant so that model A corresponds to the stationary SPDE 
model ([^, and the number of basis functions represents the number of basis functions in 
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Figure 5: Triangulation over Eastern United States; green line: U.S. boundary; red dots: 
locations of ozone monitoring stations; size proportional to the ozone levels in ppb (parts 
per billion) 


both a:-direction and ^/-direction; for example in model C, there are 3 basis functions for 
r(-) in x-direction as well as ^/-direction so that there are actually 3x3 = 9 basis functions 
for r(-). 


Table 3: Number of basis functions for the B-spline in each direction for the parameters 
K^(-) and r(-) 



A 

B 

C 

D 

E 

F 

G 

H 

I 

J 

K 

L 


1 

1 

1 

1 

1 

2 

3 

4 

5 

2 

3 

4 

^(■) 

1 

2 

3 

4 

5 

1 

1 

1 

1 

2 

3 

4 


To measure the £t and predictive performance and select the appropriate representations 
for K^{-) and t(-), we employ the leave-one-out cross validation and aim to hnd the model 
with the smallest Log Score. Figure presents the Log Scores of the 12 models for the 
two parameters k^(-) and r(-) as shown in Table [fusing BS-SPDE-G approach with d = 
1,2, 3,4, 5 respectively. It is easy to see that the Log Scores obtained from BS-SPDE-G 
with higher d are generally smaller than those obtained from BS-SPDE-G with lower d. 
Using BS-SPDE-G with a specihc d, the Log Scores for different representations of 
and r(-) are different. In general the non-stationary models B-L yield smaller Log Scores 
than the stationary model A and the overall smallest Log Score is obtained with model L 
and BS-SPDE-G d = 3. Hence model L is chosen to be the model for the non-stationarity 
in this case for further prediction. As shown in Table model L corresponds to 4 basis 
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functions for and 4 basis functions for r(-) in both x and y direction. This suggests 
that both K^(-) and r(-) display spatial variation over the domain. Note that the number 
of parameters are not taken into account for model selection here. In fact, we notice that 
the Log Score obtained using BS-SPDE-G d = 3 with model D is only slightly higher than 
model L while the number of parameters used to represent the non-stationarity is only half 
of model L. A proper model selection technique would account for it, e.g. AIC and BIG, 
but this is beyond the scope of this paper. 



Figure 6: Log Scores for the models A-L using BS-SPDE-G with d = 1,..., 5 

Then we apply the non-stationary model L for the two parameters and r and predict 
ozone levels over the Eastern United States using the BS-SPDE-G approach with d = 3. 
Figure displays the posterior mean and standard deviation of the predictions given the 
observations presented in Figure]^ As we can see, the predicted ozone level is low in the 
south-east corner and at the top of the north-east corner and high in the north and middle 
area, which is consistent with the observations. Furthermore, the posterior predictive 
standard deviation shows some spatial variation over the entire domain because of the 
irregular distribution of the observations and the non-stationarity of the SPDE model. 


4 Discussion 

We have shown that polynomial basis of order greater than one can be easily implemented 
in the SPDE framework for GFs using bivariate splines. Both the theoretical results and 
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Figure 7: (a) Posterior mean; (b) posterior standard deviation: ozone levels over Eastern 
United States predicted using BS-SPDE-G with d = 3 and model L for the two parameters 
and r(-) 


numerical simulations have demonstrated that this new approach has its advantages over 
the linear finite element approach in many applications in terms of both approximation 
accuracy and computational efficiency. By using higher degree representations, we can also 
implement the least squares solutions to the SPDE Q for a = 2. This should be more 
computationally efficient to make inference than the corresponding Galerkin solutions due 
to the different sparsity structures of the two solutions. In terms of application, we have 
shown that the SPDE approach can be applied to the spatial modelling of bathymetry. 
However, the current commonly used mapping tools, e.g. Generic Mapping Tools (GMT) 
used by NOAA (Eakins and Taylor, 2010), do not include uncertainty estimates of the maps. 


GMT also requires high smoothness conditions, see Smith and Wessel| (|1990) and Wessel 


and Bercovici (1998), that may not be appropriate for the bathymetry/topography. Hence 


the computationally efficient SPDE approach has promising potential in spatial mapping. 

There is still room for further investigation in the application of bivariate splines in the 
SPDE approach. It has been suggested by the numerical simulations that the degree of 
polynomial basis has an impact on the performance of the SPDE approach so that choosing 
appropriate degrees is essential to the SPDE approach to capture global and local spatial 
characteristics. In fact, within the framework of bivariate splines, the degree of polynomials 
can be different over different triangles 


Hu et al. 


(2007) proposed a new spline method 
which allows automatic degree raising over triangles of interest. This new method is able to 
solve linear PDEs very effectively and efficiently. Another extension to manifolds could be 
considered. Lai et al. (2009) discussed the application of spherical splines in geopotential 
approximation where the techniques of triangulated spherical splines can be applied to 
represent the Matern fields on manifolds. Furthermore, when a is larger than 2, which 
means the smoothness parameter i/ in ([^ increases as well, smoother sample paths of 
the Matern helds are expected according to [Paciorek and Schervisli (2004). However, it 
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is quite difficult to implement higher orders of smoothness in conventional finite element 
representations. But within the bivariate splines framework, higher orders of smoothness 
conditions can be implemented easily by imposing linear constraints on the B-coefficients 
(Lai and Schumaker, 2007) which leads to another potential extension of this approach 
allowing smoother representations of the GFs. However the locally supported bases for 
smoother bivariate spaces S'^(A) for r > 0 are difficult to construct and the Bayesian 
inference with large number of linear constraints are still open to research. Hence this 
needs to be investigated further in order to implement the SPDE approach with bivariate 
splines of higher smoothness. 


A Appendix 

We include some details about bivariate splines relevant to this paper. Then the sketch of 
proofs for Theorem and Propositions [^-[^ are provided. 


A.l Preliminaries on bivariate splines 

To evaluate a polynomial of degree d in B-form over any triangle, say p = J2i+j+k=d 

at the point v = {x, y) whose barycentric coordinates are b = ( 6 i, 62 , ^ 3 ) with 61 +& 2+&3 = 1 , 

let cf^l = Cijk and for all I = 1 ,..., d. 



bic 


(Z-l) 


+ 62C 


(Z-l) 

*J+1 




+ &3C 


(Z-l) 
i,j,k-\-l * 


For i + j + k = d — l,we have 


p(v) 


i+j-\-k=d—l 


for all 0 < Z < d. In particular, p(v) = CqqI- This is called the de Casteljau algorithm (Lai 


and Schumaker, 2007) 


Each vector u can be uniquely described by a triple ( 01 , 02 , 03 ) called directional co¬ 
ordinates of u, that is Oj = Qfi — /3i, i = 1,2,3, where (ai,a 2 ,a 3 ) and /32, Ps) are the 
barycentric coordinates of two points co and co such that u = lo — to. It is easy to see that 
the barycentric coordinates of a point sum to 1 , while the directional coordinates of a vector 
sum to 0. Suppose u is a vector in whose directional coordinates are a = ( 01 , 02 , 03 ), 
then for i + j + k = d, we define the directional derivative of at location v with respect 
to directional vector u to be 


C.-Bjilv) = d [oi-Stir i(v) + a 2 By-i,t(v) + 


(9) 


The integrals and inner products of the Bernstein polynomials can be calculated pre¬ 
cisely as presented in the following lemma. 
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Lemma 1 Let p = he a polynomial of degree d on triangle T (with area 

i-\-j-\-k=d 

At), then 

f At 

p{x,y)dxdy = 2^ djk- (10) 

I 2 / i+j+k=d 

Let q = Y) another polynomial of degree d on triangle T, then the inner 

U-\-^-\-K.=d 

product of p and q is 


p{x,y)q{x,y)dxdy = 


e (ii) 

\ d / \ 2 / i-\-j-\-k=d 


iy+fi-\-K,=d 

For the spline space ^^(A), the domain points are dehned to be the set 

= {^ijk = (^vi + jv2 + k^^3)/d,i+j + k = d,T = (vi,V2,V3) e A}. 
Therefore the spline function can also be denoted by 

5|t= 

5sx>d,T 


where stands for Bjf^ for ^ G V^^t and cg is the corresponding B-coefficient 

Cijk- Note that since s is continuous, if f lies on an edge shared by two different triangles 
T and T, then the corresponding coefficients cg for s\t and s\f should be the same. Then 
we show that the basis for *S'°(A) can be constructed easily with spline functions in S'^(A) 
with specific B-coefficients. For each G Pd, a, let be the spline in ^^(A) having all 
zero B-coefficients except for = 1, then we have the following result 

Lemma 2 The set of splines B = ^ forms a basis for the spline space Sj^(A) 

which satisfies > 0 and V'c(v) = 1 for all v G hi. 

It is obvious that is identically zero on all triangles that do not contain since the 
corresponding B-coefficients are all zeros so that is locally supported. 


A. 2 Proof of Theorem [T] 


The derivation of Qi 


Qf and Qc 


for a > 3 can be found directly from Appendix D.3.1 
of Lindgren et ah (2011). Here we only present briefly how to calculate the least squares 
solution and the matrix components M, K and R. 

When a = 2, by plugging the bivariate spline representation of x(u) in to the equality 
we have 


{(), - A)Tf)hWh) = (0, W), 


h=l 
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for any appropriate test functions 0. By choosing a set of test functions to be (j)h = 
— A)'iph, we have 

m 

((k^-A)' 0t,^(K^-A)ri/’sW5) = ((K^-A)V’oW^), t = (12) 

The left hand side is 


+ 2K^{V'lpt,V'lps) + {A'll!s,A'4!t))Ws, 

S = 1 


by applying the stochastic Green’s hrst identity along with the Neumann boundary condi¬ 
tions. The integral on the right hand side is in fact Gaussian with mean zero and covariance 
matrix whose {t, s)-th element is 


Gov(((k^ - A)'ijjt,W),{{K‘^ - A)^ps,W)) = + 2K^{V'4)t,V'4)s) + {A'ijjs,A'ijjt). 


Then we can write (12) in the matrix form as 

+ 2k‘^K + R)w ~ A(0, + 2k‘^K + R), 


where the (t,s)-th entry of the matrices M, K and R are respectively 
Kts = {'V4’t,'^4’s), and R^s = {A'lpt, A'lps)■ which are usually called mass matrix, stiffness 
matrix and roughness matrix in bivariate spline literature. Therefore it is easy to show 
that the precision matrix of w is Q = t^{k^M + 2k^K + R). 

Following the Lemma 0 and Vp = J2i+j+k=dCijk^Bf.^ > = T.i+j+k=dCijkABf-f^ for 

any p = J2i+j+k=d^ijk^ijk^ have the contribution of triangle T to the {t,s)-th entry of 
M, K and R for t,s = 1,..., m are 


Ktsln = V'0s)r = 

Rtsln = (V'0t, V'^s)r = c[\tRtCs\t, 

where Kt and Rt are dehned in Theorem and Ch\T is the column vector of B- 
coefficients of rjjh associated with triangle T, h = 1, Then it is followed that 

Mts = ^ Mtslr = c(MoC^, ^ Kt^lr = c(KoC^, Rt* = ^ Mtslr = c'jRqC^, 

T T T 

where Mq = diag(MT,T G A), Kq = diag(Kr,T G A) and Rq = diag(R'r,T G A). 
Therefore we have the following simple matrix representation that 

M = C'MoC, K = C'KoC, R = C'RqC. 
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A.3 Proof of Proposition [T] 

For 1 < q < oo and d > 1, the associated Sobolev space on hi in is defined by 


= ll/IU,.,n<oo}, 


where 




— / \k=0 


\f\k,q,Q 


1/g 


El/I fc,cx),n; 

k=0 


with 


, 1 < g < oo 

g = oo, 

1/9 


E mD^fWln] , 1<«<«3 

\k,q,fi — \ \u+ii=k 


max ||F)^F)(;/||oo,o, 


^ u+fi=k 


g = oo, 


and 


_ I (L l< 9 <oo, 


q,n 


esssnp^gfj |/(m)|, g = oo. 

Let /a(s) be the if ^-orthogonal projection of / G fl onto the bivariate 

spline space 5°(A), it follows that 

f f{s)LxA{s)ds = j (/(s) -/A(s))La;A(s)ds-f f f^{s)LxA{s)ds 
Jn Jn Jn 

= / /A(s)LxA(s)(is 

Jq 

= [ U{s)dW{s), 

Jn 

where the second eqnality follows from the orthogonality of /(s) — /a(s) to *S'°(A) with 
respect to inner product. Then we have 

[ fis)L{x{s) - XAis))ds= [ (/(s)-/A(s))dhF(s). 

Jn Jn 

Hence it follows from the white noise integrals that 


e( I fis)L{x{s)-XA{s))ds) =E^(/(s)-/A(s))dlT(s) 

= [ (/(s) - /A(s))^ds. 


Then it follows from standard results in bivariate splines literatures, for example Th. 5.19 
in 


Lai and Schumaker (2007) that under some suitable assumptions on the triangulation. 


we have for 1 < m < d. 


||/-/A||2,o<iL|A 


m+l I 


|?n+l 52 ,n‘ 
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A.4 Proof of Proposition 

First of all, it is easy to see that 

w'j-Mwg = (/a,5'a)a = X] / f^gAdxdy (13) 

T6A 

since /a^S'a ^ >S'°(A). Next we can see 

w'fMwg = / 4>^4>r]dxdyc^{g^) 

TeA geCd.T V 

TgA ^GT)d,T V 2 / 

where Pd,T = {(^vi + jv 2 + kw^)/d,i + j + k = d] \s the set of associated domain points 
of triangle T = (vi, V 2 , V 3 ), At is the area of triangle T and c^(s) is the B-coefficient of s. 
When /a = C is a constant C, it is easy to see that 

At 

( d-\-2 



) 


cdfA)cdgA) = / fAgAdxdy 


and hence, we have 


w}Mwg =J2f 

TeA 


fAgAdxdy 


which is Wj^luig by (13). Similar when is a piecewise constant. Also, when d = 1 , this 
resnlt follows from Lemma 1 in Chen and Thomee (1985). We now prove it for general 
d > 1. 

We hrst note that 


cdfA)cdgA) =( c 5 (/ a ) - fAiOydgA) + fAiOic^gA) -gAiO) 

+ /A(0dA(0- 


Then we claim that 




An 


id) 


fAiOgAiO approximates / f£^{x,y)g^{x,y)dxdy. 


Indeed, let ns recall the Bernstein-Bezier approximation of arbitrary continnous fnnction 
F on T. That is, nsing Th. 2.45 in Lai and Schnmaker (2007), we have 


|F-Brf(F)||T,oo < 


\Jf 

d 


2,T 


(14) 
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where Bd{F) = Bernestein-Bezier polynomials of degree d. 

Letting F{x,y) = fA{x,y)gA{x,y), we have 


F{x,y)dxdy — / Bd{F)dxdy\ 


An 


< 


IT? 

d 


fAix,y)gAix,y)dxdy - ^ /a(0^a( 0/d+2N 

5sx>d,T V 2 / 

|/Afi'A|2,rda;d|/ 


irP 

'^F —^|/A5'A|2,i,r 


<K 


d 

\Tl 

d 


|/A|2,2,T|5'A|2,2,r, 


where we have used the fact that /aS'a is a polynomial of degree 2d in the second inequality 
and the Cauchy-Schwarz inequality in the last inequality. This hnishes the proof of the 
claim. 

Next we consider 


A(r):= E 


Ai 


-(/a(0 -c^(/a))c^(£/a)- 


/ d-\-2\ 

V 2 / 

We have 

|Jl(T)| = At||{/a(0 - C^{fA)}^eVd,T\\oc\\{c^igA)}^eVa,T\\oc 
and hence, by Th. 2.6 in Lai and Schumaker| ( |2007 ), 

\hiT)\ < ATK‘^\\Bd{fA) — /aIItIs'aIt, 

where iL is a positive constant. We use the property of Bernstein-Bezier approximation 
again, i.e. the estimate in (14) to have 

itp 

\h{T)\ ^|/A|2,T||fi'A||oo,0 

<-^^||/A||2,i,r||5'A||oo,n- 

E IA(r)| < 


Therefore we have 


TeA 


Similarly, we can discuss 


l2{T)-.= E 


An 


-/a(0(cc(^a) -^a(O) 


/ d-\-2\ 

V 2 / 

to have a similar estimate as Ii{T). Putting these three estimates above we have obtained 

kA(/A,5'A)| < -^|^P(||/A||2,2,n||5'A||2,2,n + ||/a || 2,1,0 Hfi'A || oo,n + || /a || oo,0 Hfi'A || 2,l,o), 

where iP is a positive constant, |A| is the length of the longest edge in the triangulation 
A. These complete the proof. 
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